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ABSTRACT 

^ . Context. During protostellar collapse, conservation of angular momentum leads to the formation of an accretion disc. Little is known observa- 
^s^j ■ tionally about how and when the velocity field around the protostar shifts from infall-dominated to rotation-dominated. 

f--. Aims. We investigate this transition in the low-mass protostar LI 489 IRS, which is known to be embedded in a flattened, disc-like structure that 
\^ ■ shows both infall and rotation. We aim to accurately characterise the structure and composition of the envelope and its velocity field, and find 
' clues to its nature. 

1— H Methods. We construct a model for LI 489 IRS consisting of an flattened envelope and a velocity field that can vary from pure infall to pure 
rotation. We obtain best-fit parameters by comparison to 24 molecular transitions from the literature, and using a molecular excitation code 
and a Voronoi optimisation algorithm. We test the model against existing millimeter interferometric observations, near-infrared scattered light 
, imaging, and 12 CO ro-vibrational lines. 
O | Results. We find that L1489 IRS is well described by a central stellar mass of 1.3±0.4 M Q surrounded by a 0.10 M Q flattened envelope with 

6 

Conclusions. We speculate that LI 489 IRS was originally formed closer to the center of this core, but has migrated to its current position over 
• • ■ the past few times 10 5 yr, consistent with their radial velocity difference of 0.4 km s . This suggests that L1489 IRS' unusual appearance may 
. be result of its migration, and that it would appear as a 'normal' embedded protostar if it were still surrounded by an extended cloud core. 
Conversely, we hypothesize that the inner envelopes of embedded protostars resemble the rotating structure seen around L1489 IRS. 



approximate scale height h as 0.51R, inclined at 74°^\j . The velocity field is strongly dominated by rotation, with the velocity vector making an 
angle of 15° ± 6° with the azimuthal direction. Reproducing low-excitation transitions requires that the emission and absorption by the starless 
core 1' (8400 AU) east of L1489 IRS is included properly, implying that L1489 IRS is located partially behind this core. 
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1. Introduction 

The last three decades has provided a detailed understanding 
of the process of low-mass star formation through theoretical 
work and advancement s in the obser vational facilities (see for 
ex ample the revi ew by lAndre et"afl 2000; or several reviews 
in Reipu rthet alJ 2006). These achievements have given us a 
detailed view on infant stars through their various stages of 
formation. Low-mass stars form out of dark molecular clouds 
when dense regions can collapse under the influence of their 
own gravity. When sufficient density is reached, a protostellar 
object is formed, still deeply embedded in a surrounding enve- 
lope. Conservation of angular momentum leads to the forma- 
tion of a disc around the protostar onto which the surrounding 
dust and gas is accreted, although little details are known of 
how exactly discs grow. As the stellar wind starts to clear out 
the envelope, the star and the disc becomes visible in the op- 
tical and infrared and the object enters the classical T Tauri 
stage which then later evolves into a main seq uence star JShul 
ll977HLizano & Shull989tlAdams et alJll988l) . Most observed 



Young Stellar Objects (YSOs) are usually classified based on 
the shape of their spect ral energy di stribution (SED) as either 
a Class I, II, or III iLada & Wilkin J 19 84 . Class I objects are 
deeply embedded in dense cores, while Class II objects are sur- 
rounded by actively accreting discs. Class III objects have lit- 
tle material left in a disc, but are still descending to the Main 
Sequence. Sometimes, however, a YSO does not clearly fit into 
one of these categories. Those objects are most likely the ones 
that can shed light on some of the missing pieces of the picture. 
In this paper we study one such transitional object, L1489 IRS, 
and investigate the structure, dynamics, and composition of its 
circumstellar material. 

L1489 IRS (IRAS 04016+2610) is classified as a Class I 
object b ased on its SEP a nd visibility at near-infrared wave- 
lengths (Mve rs et alJll987h . Like many embedded YSOs, line 
profiles of dense gas tracers like HCO + J '-3-2 and 4-3 show 
red-shifted absorption dips usually interpreted as indications 
of inward motions in t he envelopes (Gregersen & EvansfeoOO: 
iMardones et alJ fl997l) . However, Hog erhe ijde (2001 ) shows 
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that the spatially resolved HCO + 7=1-0 emission exhibits a 
flattened, 2000 AU radius structure dominated by Keplerian 
rotation. In this aspect, L1489 IRS more closely resem bles 
a T Tauri star with a circumste l lar disc ( Koerner & Sargent 
Il995t iGuilloteau & Dutrevlll998l ISimon et al.ll2000l) . T Tauri 
discs, however, are in general much smaller than the disc struc- 
ture seen in L1489 IRS w ith radii of several h undreds of AU. 
Scattered light imaging by 1 



;ttetal] (ri999) shows the cen- 
tral stellar object and the presence of a slightly flaring dark 
lane, consistent with the disc-like configuration inferred from 
the HCO + 1-0 obse rvations. Care ful analysis if the circumstel- 
lar velocity field bvlHo gerheiidel {2001 ) revealed that infalling 
motions are present at ~ 10% of the Keplerian velocities. 
iHoeerheii de hypothesized that L1489 IRS is in a short-lived 
transitional stage between a collapsing envelope (Class I) and 
a viscously supported, Keplerian disc (Class II). Observations 
of ro-vibrational CO absorption lines at 4.1/j.m showed that the 
inward motions cont inue to within 1 AU from the central star 
feoogert et al.l2002h . 



Table 1. Single-dish Molecular Line Data Set 



In this paper we address a number of questions about 
L1489 IRS. We construct a model for the circumstellar struc- 
ture that accommodates all observations, ranging from an ex- 
tensive set of single-dish molecular line measurements to the 
interferometric observations, the scattered ligh t imagi n g, and 
the CO ro-vibrational absorption lines. Hoger heiidd d200ll) 
adopted a flared disc model with a fixed scale height for the 
structure inspired by the interferometric imaging. In this Paper 
we choose a description for the circumstellar structure that 
can be smoothly varied from spherical to highly flattened, 
and investigate if the full data set requires a disc-like config- 
uration. We also adopt a velocity field that can range from 
purely Keplerian to completely free-fall, or any combination of 
the two. By considering the full data set, stronger constraints 
can be set on the velocity field and the dynamical state of 
L1489 IRS than possible before. We perform a rigorous opti- 
misation of the model for L1489 IRS using all available single- 
dish line data, and test the model by comparing the interfer- 
ometric observations, the scattered light imaging, and the CO 
ro-vibrational absorption lines to predictions from the model. 
Once we have established a satisfactory model, also taking into 
account the immediate cloud environment, we explore the na- 
ture of L1489 IRS. Does it represent a transitional state be- 
tween Class I and II? Do all YSOs go through this stage? Or is 
L1489 IRS is some way special? 



The layout of this Paper is as follows. Section |2] present 
our data set and a detailed overview of the model and fit op- 
timisation procedure. Section [3] describes our best fit model 
and its reliability, including comparison with observations of 
interferometric observations, scattered light imaging, and CO 
ro-vibrational absorption lines. Section0]discusses our results 
in the light of the nature and evolutionary state of L489 IRS 
and explores the wider implications for our understanding of 
star formation. Section|5]concludes the Paper with a brief sum- 
mary. 
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1.5 kms- 1 . 

The error bars on the intensity does not include the 20% calibration 
error. 



2. Observations and model description 

2.1. Single-dish observations 

The prim ary data set on L14 89 IRS use d in this paper was pub - 
lished by Hog erheiide et al.ldl997l) and Ltorgensen et al.l lEbo4). 
and consists of 24 transitions among 12 molecular species. 
Figure[0shows all 24 spectra. Table^lists the transitions, in- 
tegrated line strengths, line widths, and relevant beam sizes of 
the single-dish telescopes. In all cases, line intensities are on 
the main-beam antenna temperature scale, using the appropri- 
ate beam efficiencies. The integrated intensities are obtained 
by fitting a Gaussian to the line. In some cases, no lines are 
visible above the noise level, and 3cr upper limits are given. 
The signal-to-noise ratio of the HNC 4-3 and H 2 CO 5i 5 -4i 4 
spectra was insufficient for a proper Gaussian fit; instead the 
spectra are simply integrated from -4 to +4 km s with re- 
spect to the systemic velocity of +7.2 km s _1 . In addition to 
these molecular line data, we also use the total mass derived 
from the 850 um conti nuum observations by JCMT/SCUBA 
(Hoge rheiide & S andelll2000h . 

Apart from the single dish data which we use in the model 
optimisation, we compare predictions by the model to other 
previously published observations: a HCO + 7=1-0 interferom- 
eter m ap from the BIMA and OVRO arrays (Hoger heiide et al] 
1998), CO ro-vibrational absorption spectra from the Keck 
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Fig. 1. The 24 single-dish molecular line spectra used for the model optimisation (histograms). The solid lines show the best-fit 
results for the model of L1489 IRS (See also Fig.EJl. 
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Telescope! Bo ogert et al.l2002l) . and the near- infrared sc attered 
light imaging from HST/NICMOS JPadeett et all 19991) . 

2.2. The model 

Previous studies of L1489 IRS clearly indicate that an axi- 
symmetric description of its circumstellar structure is required. 
In the following subsections we construct a description of 
the density n(r, 0), gas temperature T(r, 0), and velocity field 
v(r, #)=(v's, v z , vty). Throughout we attempt to keep the number 
of free parameters at a minimum. In the end we arrive at four 
free parameters, in addition to the eight molecular abundances 
which we fit but assume constant throughout the source, and 
seven parameters that we hold fixed (Table [2}. 

2.2.1. Density 

We adopt an axi-symmetric description of the gas den- 
sity n(R, z) consiste nt with the spherical model from 
( JOr gensen et alJEoO^ . These authors deduce a total mass of 
0.097 M and a density following a radial power Law with 
slope -1 .8 between radii of 7.8 and 9360 AU. We truncate this 
model at the observed outer radius of L1489 IRS of 2000 AU, 
but keep the power-law slope and mass conserved. Instead of 
a simple radial power law, n oc r~ p with p = 1.8, we adopt 
a Plummer-like profile, n oc [1 + (r/ro) 2 ] - ^ 2 with ro=4.0 AU. 
This description keeps the density finite at all radii, but since 
ro is much smaller than the scales of interest h ere, the resu lting 
densi ty distribution is identical to that used by (Jorg ensen el alJ 
2002). 

From this spherically symmetric density distribution we 
construct an axi-symmetric, flattened configuration by multi- 
plying by a factor sin^ 0, where / can take any value > (see 
Stamatel los et alJ2 004. where this approach was used for mod- 
elling starless cores). The adopted density distribution now be- 
comes, 

n(r, 0) = njyl + (-^) 2 ) ^ sir/ 0. (1) 

For / = this reduces to a spherically symmetric struc- 
ture, while for / > 10 the resulting profiles becomes largely 
indistinguishable as they approach a step function (Fig.|2}- The 
mass contained in the structure is kept constant at 0.097 M 
by adjusting no as / is varied. The only free parameter in the 
density description is the flattening parameter /. 

2.2.2. Temperature 

The temperature of the gas and the dust (which we assume to 
be identical) in the circumstellar structure of L I 489 IRS de- 
pends on the stellar luminosity which is ~ 3.7 L Q dKenvon et alJ 
Il993al) and the infrared radiative transfer through the structure. 
Since most of the circumstellar material is optically thin to far- 
infrared radiation, the deviations introduced by the flattening 
on the temperature structure are minor. Furthermore, the line 
excitation does not depend strongly on small temperature dif- 
ferences. A spherically symmetric description of the tempera- 
ture therefore suffices. Using the continuum radiation transfer 
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Fig. 2. Progression of flattening of the adopted density structure 
as / is increased from (purely spherical) to 10 in Eq. 0. 

code DUSTY (Nen kova et al.lll999l) and the density structure 
of Eq. with p = 1.8 and / = 0, we find that the temperature 
is well described by, 

/ r \-°- 35 

T(r) = 19.42 K( . (2) 

U000AU/ 

There are no free parameters in this description of the tem- 
perature. 

2.2.3. Velocity field 

The velocity field is parameterized by a central, stellar mass 
M+ and an angle a in such a way that, 

vr =-V2 ^/^sina, (3) 
V0 = ^fcK cos a. (4) 

For a — this reduces to pure Keplerian rotation around a 
mass without any inward motions; for a — | the velocity 
field is that of free fall to a mass M*. Intermediate values of 
a produce a velocity field where material spirals inward. The 
implicit assumption in this description is that both components 
of the velocity field vary inversely proportional with yfr. Note 
that a should not be confused with the geometric angle deter- 
mining the direction of the flow lines. 

In this description there are two free parameters, the stellar 
mass M+ and the angle a which is kept constant with radius. 
In addition to this ordered velocity field, we add a turbulent 
velocity field with FWHM 0.2 km s" 1 . 

2.3. Molecular excitation and line radiative transfer 

The excitation of the molecules and the line radiative transfer is 
calculated using the Accelerated M onte Carlo code RATRAN 
(Hogerheiide & van derTakl2 000). Collisional excitation rates 
are taken from the Leiden A tomic and Molecular Database 
LAMDA ( Sch oier et al.l2005h . We lay out the model onto three 
nested 8x6 grids (Fig.|3j- The innermost grid cell is subdivided 
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Fig. 3. Layout of the grid cells for a model with / = 3.8. 



four times, so that the innermost cell is resolved down to 4 AU. 
All properties are calculated as cell averages, by numerically 
integration over the cell and divide by its volume. To reduce 
computing time, cells with H2 densities below 10 3 cirT 3 are 
dropped. Such cells do not contribute significantly to the line 
emission or absorption. Dust continuum emission is included 
through a standar d gas-to-dust rat i o of 1 00:1 and dust emis- 
sivity from Ossenkotrf & Henning ( 1994) for thin ice mantles 
which has been accreted and coagulated for about 10 5 years. 

Synthetic observations are created from the molecular ex- 
citation by performing ray-tracing after placing the object at a 
distance of 140 pc and an inclination i (a free parameter). The 
resulting spectra are convolved with the appropriate Gaussian 
beams. Figure^shows the best-fit model spectra (the best fit is 
discussed in section[3}. 

2.4. Modelling the neighbouring cloud core 

During the optimisation of the fit ( 32.51 it became obvious 
that several lines, and especially those of lower-lying rotational 
transitions taken in large beams, were contaminated by emis- 
sion with small line width. This emission component is espe- 
cially clear in the C 18 1-0 and 2-1 lines, the CS 2-1 line, the 
HCN 1-0 line, and, to some extent, the HCO + 1-0 line (Fig.Q). 
The emission has a Vlsr of 6.8 km s _1 , slightly lower than that 
of L1489 IRS of 7.2 km s _1 . Cold fore- or background gas with 
small turbulent velocity is the li kely cause for this component . 
The 850 fim SCUBA map from lHogerheiide & Sandeif lEoOO) 
reveals that L1489 IRS sits at the edge of an extended, probably 
starless, cloud core with a radius of 60" (8400 AU). Cold gas 
in this core therefore contributes to the low-/ emission lines, 
and especially in spectra taken with large beams. 

We construct a simple description for the neighbouring 
core, so that we can take its emission into account in our opti- 
misation of the model for L1489 IRS, as well as its absorption 
if this source is located behind the core. We approximate the 
core as spherical with a radius of 60", which is roughly the dis- 
tance of L1489 IRS to its centre. We assume that it is isother- 
mal at 10 K and that it has abundances typical for starless cores 



( Jeir gensen et al.lEo04h . For the species which does not show 
any cloud core emission, the abundances are unconstrained and 
we just set the abundances sufficiently low. In the case of CO 
we use an abundance of 5 x 10~ 5 . The CS abundance is set to 
2 x 10~ 9 , and the HCO + and HCN abundances are 27 x 10~ 9 and 
4 x 10~ 9 respectively. We deri ve its density distribution by fit - 
ting the 850 fim emission from Hogerheiide & Sandell (2000). 
We find an adequate fit for a radial power-law with slope -2 
and a density of 4 x 10 6 crrT 3 at r = 1000 AU resulting in a 
cloud mass of 2.9 M . This is consistent with the drop off in 
density found in many starless cores on scales (r > 1000 AU) 
that are relevant to us dAndre et alJl996|) . Because it falls out- 
side even our largest beam on L1489 IRS we do not investigate 
if the density in the neighbouring core levels off at the center, 
as is seen for many starless cores. The relative smoothness of 
the 850 fim emission suggest that this is the case, however. 

Using RATRAN we calculate the expected emission and 
the optical depth of each of the observed transitions. In our 
model optimisation procedure (see below), the emission from 
L1489 IRS and the neighbouring core are added on a channel- 
by-channel basis, with the appropriate spatial offset for the 
core. We find that we can only make a fit that is reasonable if 
LI 489 IRS is located behind the core; we need both the emis- 
sion and the opacity of the cloud. This is taken into account 
by first attenuating the emission from L1489 IRS by the core's 
opacity, again on a channel -by-channel basis, and subsequently 
adding the core's emission in each channel, followed by beam 
convolution. 

In this section we derived only an approximate model for 
the neighbouring core. Its effects are taken into account in the 
model spectra, but the description of the core is not accurate 
enough to include in the model optimisation. This would re- 
quire a much more detailed analysis than possible here. In the 
procedure outlined in the next section, we therefore mask out 
those regions in the spectra strongly affected by the emission 
and absorption of the core. 

2.5. Optimising the fit 

Our model has four free parameters: the inclination i, the flat- 
ness parameter /, the stellar mass M*, and the angle of the 
velocity field a. In addition, the abundances of the molecules 
are unknown. All other parameters are held fixed. Table|5]lists 
the parameters. 

Considering the size of the parameter space and the time it 
takes to calculate a single spectrum 1 the task of finding the pa- 
rameter vector resulting in the best fit is non-trivial. This is fur- 
ther complicated by the degeneracy of the model results to dif- 
ferent parameters. For example, increasing the abundance can 
have the same effect on the line intensity as increasing the incli- 
nation or the flatness, but these will have very different effects 
on the line profile shape. 

Instead of calculating all possible models in the allowed 
parameter space, we use Voronoi tessellation of the parameter 

1 Depending on the species and the optical thickness, we can cal- 
culate a spectrum in between five minutes and half of an hour, on a 
state-of-the-art desktop processor. 
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cube (see e.g. lKiangll966l for details on Voronoi tessellation). 
A random set of n points p„ in the parameter cube is picked and 
model spectra are calculated for each of these. Then the param- 
eter cube is divided into Voronoi cells, defined as the volume 
around a point in the parameter cube containing all points 
q closer to pi than to any other of the points p n (n + i). The 
parameter cube is scaled in arbitrary units, so that the allowed 
parameter ranges falls between and 1 . On this dimensionless 
unit cube a simple metric in d dimensions is used to define the 
cells, 



S^^q.-Pif 



(5) 



assuming that the solution depend linearly on all parameters. 
This assumption is not true especially for large values of s, but 
because we have no knowledge of the geometry of the param- 
eter space, we use the simplest possible measure. In order to 
minimise the effect this has on our final solution we can in- 
crease the initial sample rate so that the average distance be- 
tween the points becomes smaller. After one or two iterations, 
the volume of each cell is small enough so that the assump- 
tion of linear dependence is good. By scaling the parameters to 
the same range we make sure that each parameter is weighted 
equally in the distance measure. 

The cell which contains the point pi resulting in the best fit 
is chosen, and a new set of random points are picked within this 
cell, and the procedure is iterated until sufficient convergence 
has been achieved. This method is only guaranteed to reach the 
true best fit if only one global minimum exist and if there are 
no (or few) local minima. To check whether we find the true 
optimum, we make several runs, with different randomly dis- 
tribution initial points. We find that we always reach the same 
minimum, and conclude that local minima are few and not very 
deep. 

For every calculated model spectrum, the fitness is evalu- 
ated by regriding the model spectrum to the channel width of 
the corresponding observed spectrum, centering it on the LSR 
velocity of 7.2 kms~' , and calculating the^ 2 between the model 
and the observed spectrum, 



= -7-7 

M ^ Nm ^ 



(I(ri)obs - I(n) mo del) 2 



(6) 



where M is the number of spectra and N is the number of veloc- 
ity channels in the m'th spectrum. This way we give an equal 
weight to all spectra even though the number of channels vary 
in each spectrum. Those channels affected by the neighbouring 
core are not included in the x 2 measure. Every spectra has a 
fixed passband of 14 kms~' so that an equal amount of baseline 
is included for each spectrum. 

Using this method, with a set of 24 random points per it- 
eration, we converge on an optimal solution after four to five 
iterations, corresponding to 10 to 12 days of CPU time. For 
practical reasons we initially chose only to consider the most 
structured lines (CO, HCO + and CS), lowering the computa- 
tional time to about a single day and getting a quick but rough 
handle on the initial parameter cube. We then included the other 
lines to obtain the overall best solution. 
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Fig. 4. The histograms show the x 2 distributions of models 
around the best fit position where only one parameter is varied 
at each time. A Gaussian, centered on the best fit in each panel, 
is fitted to the distributions. The dispersion of the gaussians is 
given in each panel. 

2.6. Error estimates 

Getting a handle on the uncertainties in the obtained parameter 
values is a difficult matter due to the size and complexity of the 
parameter space. As mentioned above, we have no knowledge 
of the overall geometry of the parameter space and given the 
long computation time of the optimisation algorithm, we can- 
not make a correlation analysis of each pair of parameter and 
neither can we make x 2 surfaces. Still, it is very important to 
get an estimate on the stability and reliability of our solution. 

A simple error analysis is done for the four model de- 
pendent parameters, the flatness, the velocity angle, the stellar 
mass, and the inclination, by fixing three of the parameters at 
their best fit values, and calculating models in which the fourth 
parameter is gradually increased from its lower boundary to 
the upper boundary. Histograms of the resulting (inverse) x 2 
values is shown in Fig.0] The^ 2 values are approximately nor- 
mally distributed, with the main discrepancy in the high values 
of the inclination, velocity field, and the flatness. This relates to 
the non-linear nature of the trigonometric functions associated 
with these three parameters. 

A Gaussian has been fitted to each of the histograms in 
FiglU The centre of the Gaussian is fixed on the best fit value 
and the hight is fixed by the x 2 value of the best fit so that only 
the variance, cr 2 , is free. Reasonable fits are achieved for each 
parameter with the cr value given in each panel. These values 
are taken to be a rough estimate of the magnitude of the error 
in each of the four parameters. For the inclination and flatness, 
where the error is greater than the allowed parameter range, 
the error is of course determined by the physical constrains on 
the parameter value (e.g., the inclination cannot be greater than 
90°). Note that the error bars are typically smaller than the ex- 
plored range in each parameter by a factor of 2-10. 
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With this kind of one dimensional error analysis we do not 
take into account the fact that the parameters are likely to be 
highly correlated. A few exploratory calculations, where one 
parameter was held fixed at its best fit value while the other 
three where randomly pertubed around their best fit values, 
indicated that there indeed exist a strong correlation between 
the parameters. Indeed, a degeneracy exists between the cen- 
tral mass and the inclination, which again is degenerate with 
the flattening. Only a full parameter space study can fully dis- 
entangle this and is beyond the scope of this paper. 

Because the abundance parameter mainly serves to scale 
the intensity in every channel of the spectrum, and does not 
change the shape of the profile much, this kind of error analysis 
is of little use. For any combination of the four free parameters 
that reproduces the observations, a corresponding abundance is 
found from the optically thin isotopic lines. These abundances 
are relatively insensitive to the exact geometry because of their 
optically thin nature. Therefore we assume that the error in the 
abundance values obtained here is entirely dominated by the 
20% calibration error of the observed spectra. 

Throughout this work we have assumed constant abun- 
dance for all the molecular species. In reality, abundances will 
depend on the chemistry and molecules will freeze out below a 
certain temperature. This gives rise to a drop in the abundances 
at a certain radius and it will affect, to some extent, the shape of 
the profiles but more prominently, the line ratios. Specifically, 
by removing low temperature material from the gas phase, low 
excitation lines become relatively weaker. Our model does not 
suffer from the problem of over-producing the low J lines, ex- 
cept for the case of HCO + ; a more complex abundance model 
would likely provide a better fit to the J - 1-0 and 3-2 lines. 
However, this would require a careful chemical analysis which 
is beyond the scope of this paper. A few tests showed that let- 
ting CO freeze out at 20 K does not change the best fit param- 
eters significantly, except for the abundance which will then 
have to be re-optimised. 

3. Results 

Figure \l\ compares the data to the synthetic spectra based on 
the best fit model obtained with the optimisation procedure 
described above. The results for the combined emission of 
L1489 IRS and the neighbour core is shown in Fig. [5] Because 
the emission from the core only contributes to the low J lines, 
this figure only shows the species in which the combined spec- 
trum show any difference from the L1489 IRS spectrum alone. 
For all species not shown in Fig. [5] the combined spectra is in- 
distinguishable from the one shown in Fig.^ Table[2]lists the 
parameters of the best fit. 

The inclination angle of 74° falls within the range of 60° 
to 90° which is in ferred from the scattered light image of 
Pad gett et alJ i 19991) and the modelling o f the infrared spectral 
energy distribution iKenvon et alJ dl993bl) . Section l3.2.2l shows 
that this inclination and the flattening parameter / = 3.8 re- 
produce the scattered light image, including the detectability 
of the central star. The resulting density distribution can be 
well approximated by a disk with a vertical density distribu- 
tion oc e - z ~l 2h ~ and a scale height of h m 0.57/?. The maxi- 



Table 2. Best Fit Parameters 



Free Parameter 


Value 


Inclination of disc, i 




Flatness of disc, sin^(6') 


J.O 


Central mass, M, 


1.35 Mq 


Velocity angle, a 


15° 


Abundance of CO" 


3.5 x 10~ 5 


Abundance of CS 


0.5 x 10- 9 


Abundance of SO 


2.0 x 10-* 


Abundance of HCO + 


1.9 x 10 -9 


Abundance of HCN 


0.2 x 10~ 9 


Abundance of HNC 


0.2 x 10" 9 


Abundance of CN 


0.2 x 10-" 


Abundance of H 2 CO 


0.7 x 10-" 


Fixed parameters 


Disc mass 


0.097 M Q 


Temperature slope 


-0.35 


Temperature at 1000 AU 


19.42 K 


Density slope 


-1.8 


Turbulent velocity, FWHM 


0.2 km s" 1 


Outer radius 


2000 AU 


Distance 


140 pc 



" The main isotopic abundance has been derive d from the C ls O abun- 
dance using a 16 0/ 18 ratio of 540 JWilson & RooJl994 . 



mum deviation of this approxim ation is only 4%, up to an an- 
gle of 60° above the midplane. iHogerheiida J200ll) described 
L1489 IRS with flared disc with an adopted density scale height 
of h — 0.5/?, so our best fit model shows that the structure 
mig ht in fact be fl atter than previously assumed. However, 
Hogerheijde (2001) used an inclination of 90° for the flared 
disc (a flat disc at i = 60° was also tried), so the projected 
column density distribution are quite similar in both cases. 

The best fitting velocity vector makes an angle of only 
15° with respect to the azimuthal direction, consistent with 
Hog erheiidd (fcoOl) who shows that rotation is the primary 
component in the velocity field. However, our central stellar 
mass of 1.35 M n is c onsiderably higher than the 0.65 M Q de- 
rived bv lHogerheiidej Part of this is due to different definitions 
of the velocity field. iHogerheiide reports a mass expected for 
Keplerian motion based on the azimuthal component of the ve- 
locity field alone, while the mass derived here results in cos a 
times the Keplerian velocity (Eq. 0}. In addition, a different 
inclination is found. These two factors together would give a 
mass of 1.2 M p for our mod el. This mass is still 80% higher 
than that from Hogerheijde but our best fit has a y 1 value of 
1.69. That is nearly half that of the best model ofllHo gerheiidel 
of ~3. We ascribe this difference to a more thorough search of 
the parameter space. 

The best fitting abundances are consistent with the a bun- 
dances obtained for L1489 IRS in (Jeirge nsen et alJl2004l) to 
within a factor of 2-3 although the values obtained here are all 
higher. This may be due to higher line opacities i n our model as 
a result of the different adopted velocity fields; J0rg ensen et alJ 
do not include a systematic velocity field and only a single tur- 
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Fig. 5. The combined emission from the L1489 IRS model and the neighbour cloud core. Only spectra which are affected signif- 
icantly by the cloud core emission are shown here. 



bulent line width. Consistent with this previous work, we find 
no evidence for depletion of CO, in accord with the relatively 
high temperatures exceeding the 20 K evaporation temperature 
of CO throughout most of the disc. Another interesting finding 
is the CN/HCN abundance ratio of 1 .0, which is more reminis- 
cent of dark clouds than of circumstellar discs, suggesting that 
chemically, L1489 IRS is close to its original cloud core and 
that p hoto-dissociation does not play a major role yet jThi et alJ 
2004). 



3.1. Quality of the best fit 

The overall correspondence of our model to the data is good 
and most of the spectra are well reproduced. The line widths 
of C 18 and C n O are well reproduced, and the C n O 2-1 and 
3-2 lines are found to exclusively trace L1489 IRS and be un- 
contaminated by the neighbouring core. All three C 18 and 
the C 17 1-0 line have narrow line peaks that originate in the 
neighbour core; some C I8 J = 1-0 emission that is not repro- 
duced can either be caused by additional material along the line 
of sight, or be due to the approximate nature of our description 
of the core. 

The sulfur bearing lines are not very intense and show little 
structure. Again we see a narrow peak in the CS J = 2-1 al- 
most entirely accounted for by the neighbouring core and per- 
haps also in the SO line. The non-detection of C 34 S 2-1 places 
an upper limit on the CS abundance in the neighbouring cloud 
of 2 x 10~ 9 . The CS 7-6 line is poorly fit, but the observed 
spectrum has low signal-to-noise. 

The HCO + lines are the strongest among our sample and 
show most structure in their profiles. It was important to be 
able to reproduce the double peak in the HCO + J = 4-3 line 
because this feature is a very clear tracer of the velocity struc- 



ture in L1489 IRS. Our model is able to reproduce this feature. 
However, we do not observe the double peak in the HCO + J - 
3-2 line but only a slight asymmetry. This provided a major 
constraint on the velocity field. The neighbour core does not 
contribute in the 4-3 line but it makes up for almost all of the 
excess emission seen in the 1-0 and 3-2 lines of HCO + . 

Of the nitrogen bearing species, HCN 7 = 1-0 with its three 
hyperfine components shows very narrow lines, which is well 
reproduced by the neighbour core model which dominates the 
emission. The HCN J = 4-3 line is much broader and uncon- 
taminated by the neighbouring core. A narrow peak in the spec- 
trum at 8 km s cannot be due to the neighbour core and we 
assume it is noise. The HNC 4-3 line has a rather low signal 
to noise ratio and can only provide a reliable estimate within a 
factor of a few. CN J = 1-0 also has a low signal to noise ratio, 
but in the case of CN, the abundance is well constrained by the 
J = 3-2 line. 

All spectra discussed above have been obtained with single- 
dish telescopes that do not resolve the 2000 AU radius source. 
One can wonder if it is justifiable to use a non-spherical model 
when the single-dish observations does not contain any spa- 
tial information. Although the scattered-light and interferome- 
ter images show that the structure is non-spherical, one could 
argue that introducing the flattening is simply a way to improve 
the fit by adding a free parameter. However, this is not the case. 
Figure [6] shows the best fit to the HCO + J = 4-3 line, with a 
spherical model (/ = 0) and a model where the flattening and 
inclination are free parameters. The flattened model provides 
a considerably better fit with lower x 1 value. And, more im- 
portantly, the spherical model shows too much self absorption. 
Our optimisation algorithm returns the model which minimises 
the difference between data and model for each velocity chan- 
nel. In a spherical model the column simply becomes too high. 
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Fig. 6. a) Best possible fit with a spherical model to the HCO + 
J = 4-3 line, b) The best fit to the same line where the flattening 
/ and the inclination i are kept as free parameters. Below, the 
residuals are shown with the mean and two standard deviations 
indicated. 



By flattening the model and adjusting the inclination, we can 
exactly reproduce the right amount of self absorption seen in 
the data while keeping the total column density high enough to 
produce the right line strength. The degeneracy between these 
two parameters are resolved by the velocity field, and thus it 
turns out that it is actually possible to retrieve spatial informa- 

tio n from single dish observation s. 

IWard -Thompson & Bucklevl J200ll) have argued that the 
amount of self-absorption can also be regulated by adjusting 
the turbulent velocity dispersion. We can indeed change the 
quality of the spherical fit by changing the amount of turbu- 
lence. We cannot, however, do that without also changing the 
velocity distance between the two peaks. Thus in order to fit 
the line width we must decrease the mass and thereby the mag- 
nitude of the velocity field, which no longer reproduces the ob- 
served infall asymmetry. We therefore find that varying the tur- 
bulent velocity width in a spherical model does not reproduce 
the observations. 

3.2. Comparison to other observations 

With the fit parameters derived above, we can now test our 
model by comparison to other observations of L1489 IRS not 
used in the fit. The neighbouring cloud is not considered in the 
following. 

3.2.1. Interferometer image of HCO+ 1-0 

We have used our model to produce a synthetic interferometer 
map of the HCO + 7=1-0 emissio n which can be dir ectly com- 
pared to the data presented by iHogerheiidd (1200 II) . Figure 
compares the model predictions of the integrated intensity and 
velocity centroid maps to the observations reproduced from 




Fig. 7. a) Interferometer map in HCO + J = 1-0. b) The corre- 
sponding synthesised map based on our best model, c) and d) 
Models with / = 1 and / = 8 respectively. The white contour 
lines show the integrated intensity, starting at 0.25 Jy birr 1 and 
increasing in steps of 0.5 Jy birr 1 . The colour scale shows the 
velocity centroid in units of km s _1 . 



Hogerheijdi feOOll) . We also show a model with / = 1 and 
on with / = 8 for comparison. The model images were made 
by taking the unconvolved image cubes from RATRAN and 
then, using the (u, v) settings from the original data set, making 
synthetic visibilities with the 'uvmodel' task from the MIRIAD 
software package. This is what we would get if the interferom- 
eter observed our model object. Then we applied the usual de- 
convolution with the invert, clean, and restore routines in order 
to reconstruct an image from the visibilities. 

The resulting synthetic image of our best fit model resem- 
bles the observations closely within the uncertainty of the abun- 
dance which sensitively affects the apparent size, when it is 
taken into account that the observations also partially recover 
the neighbouring core that is ignored in the model (Fig.0 panel 
b). However, the / = 8 model in panel d) is also in good agree- 
ment with the data which partially can be explained by the fact 
that, due to the non-linearity of the sine function, the difference 
between an / = 3 and / = 10 model is much less than the 
difference between an / = 1 and / — 3 model. This is also 
reflected in Fig |2] In any case, panel c) is obviously in poor 
agreement with the data, which shows that in order to fit the 
interferometer data, a flattened structure is needed. 

This kind of analysis is very useful to investigate the spa- 
tial distribution of the emission which is lacking in the single 
dish data. Importantly, we see that the flattening which we in- 
troduced and the amount of which we determined from the line 
profile results in a projected shape that is very close to what we 
see in the data image. Also, because the cuts are the same in 
both panels, the extent of the emission and therefore the physi- 
cal size scale of the model is consistent. 
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Fig. 8. The top left pane l shows the false co lour NICMOS im- 
age of L1489 IRS from lPadgettetalJ Jl999t) . In the three re- 
maining panels are shown synthetic three-colour composite im- 
ages based on our model. 



3.2.2. Near-infrared image 

We subsequently test our best-fit model through a comparison 
with t he near-infrared scattered light image of IPadgett et alJ 
Jl999l) . A reproduction of this image is shown in the top left 
panel of Fig. [8] Using the density structure and inclination 
found in section [5] and the properties of the central star as de- 
scribed in the introduction as input, we calculated the scatter- 
ing light emis sion with RADMC, a two-dim ensional radiation 
transfer code JPullemond & pominikll2004l). Then, using the 
ray-tracing code RADICAL JPullemond & Turollal EoOQ). we 
extract the fluxes over the full spectral range of 1 to 850 //m and 
image our model in the F110W, F160W, F210W NICMOS 
bands. 

Our best fit model of section [3] automatically provides a 
good fit to the sub-millimetre part of the spectral energy dis- 
tribution, confirmi ng that our model is consistent with the re- 
sults of Jcirge nsen et al.U2002l) . The near-infrared observations 
are more difficult to match. The column density in the inner 
part is too high to provide the clear view of the central star that 
shows prominently in NICMOS images. However, if we reduce 
the scale height in the inner 250 AU to h — 0.15/?, the col- 
umn density is reduced and the central star becomes detectable 
in the near-infrared. This adaptation does not affect the sub- 
millimetre emission or the molecular line intensities viewed 
in the much larger single-dish beams. Recently obtained high- 
resolution data from the Submillimeter Array (SMA) in Hawaii 
may be able to test this assumption of the geometry in the inner 
2"= 280 AU (Brinch et al. in prep.). 

The resulting three-colour composite images of our best fit 
model as well as a f=l and f=8 model, is shown in Fig. [S] 
Although this can only be a qualitative comparison, our model 



is able to reproduce most of the striking features evidenced by 
the observations. The opening angle of the dust cavity is found 
to be somewhat dependent on the flattening parameter / in our 
model. Although it is difficult to judge the agreement, this re- 
sult seems to favour a relatively low value of / as opposed to 
what we found above for the interferometric map. The non- 
monotonic behaviour along the axis in the / = 1 model is a 
combination of the finite gridding of the density structure and 
the effect of scattering. It has a very narrow cavity and so the 
base of the scattering nebula is actually absorbed by the high 
density material in the inner part. In the other models the cavity 
is large enough to allow all the scattered photons to escape. 

We also reproduce the near-infrared colours, confirming 
our adopted density distribution between ~100 AU and ~2000 
AU. Each of the NICMOS fluxes are reproduces to within about 
40%. Although a detailed modelling of L1489 IRS on these 
scales where most of the near-infrared emission is coming from 
was beyond the scope of this paper, we present this prediction 
to demonstrate that it is possible to combine the information 
contained in near-infrared images and sub-millimetre single 
dish measurements to obtain a self-consistent model on scales 
ranging from within a few AU up to several thousand AU. 

The cavity seen in the NICMOS imag e is likely associated 
with a molecular outflow. However, Hogerheii de et alJ {l998) 
show that L1489 IRS only drives a modest 12 CO outflow, and 
little or no impact on the line profiles is expected. Therefore, 
an outflow has not been incorporated into the model described 
in this paper. 



3.2.3. 12 CO 4.7/im absorption bands 

Finally we apply our model to fit the CO ro-vibrational ab- 
sorption bands, again by using RATRAN. Here we assume that 
initially all CO molecules are in the vibrational ground state 
(v=0) but can be excited to a v=l state by absorption of a 
photon from the central star, producing the absorption lines in 
the P and R branches corresponding to AJ = +1 transitions. 
Observations of these bands from the Keck/NIRSPEC instru- 
ment has been presented bv iBoogert et alJ J2002I) . The spec- 
trally resolved absorption lines reveal ed inward moti ons up to 
100 km s -1 . Using the same method as lBoogert iHDwecalcu- 
lated the ro-vibrational absorption lines in our model, and plot 
the average of the P(6)-P(15) lines in Fig. [9] 

We find that ou r model fits the data as well as the results 
from Boogert et all who use a contracting flared disc with 
power laws for the temperature, density, and infall velocity. 
This means that the amount of inward velocity that we ob- 
tained from the fit to the single dish lines together with the 
adopted density profile can explain the observed infall, even 
at radii much smaller than probed with the single dish lines. In 
Fig.[9]we see that absorption is present all the way up to at least 
100 km s _1 . In our model, this velocity translates into a radius 
of 0.02 AU. Material absorbing at 50 and 20 km s _1 is located 
at 0.06 and 0.4 AU respectively. 
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Fig. 9. An average of the observed 12 CO P(6)-P(15) lines with 
a similarly averaged model result over plotted. 



4. Discussion 

In this paper we derive an accurate model for the circumstellar 
material of L1489 IRS. We find that it is well described by a 
flattened structure with a radius of 2000 AU, in sub-Keplerian 
motion around a 1.35 M Q central star. While t he structure re- 
sembles discs found around T Tauri stars JSimon et alJl200ol 
e.g.), its 2000 AU radius is much larger than T Tauri discs (typ- 
ically several hundred AU). Also, discs around T Tauri stars 
are often well described by pure Keplerian motions except for 
a few cases, e.g. AB Aur, in which outwards non- Keplerian 
motions are measured (Piet u et alJEool Lin et alJkOOd) . In 
this section we discuss the evolutionary state of L1489 IRS, 
and particularly whether it represents a unique case or if other 
forming stars may go through a similar stage. We start by deriv- 
ing the life span of the current configuration. We then explore 
the relation of L1489 IRS to its neighbouring core. Finally, we 
discuss a number of open questions that only new observations 
can answer. 

By integrating the trajectory of a particle at 2000 AU we 
find the infall time scale to be 2.3 x 10 4 years. Dividing the 
0.097 M of the circumstellar material yields a mass accretion 
rate of 4.3 x 10~ 6 M yr -1 . Estimating the radius of the central 
star from the mass-radius relation (Stahler Hl988l) to be 4 R Q 
results on an accretion luminosity of L [ICC =GM,M/R, * 46 L 
corresponding to about t en times the observe d bolometric lumi- 
nosity from the YSO lKenvon et all Jl993ah . This suggest that 
not all inspiraling material falls directly ont o the star. The resul t 
is somewhat higher than the one found by Hogerheiide ( 2001 ), 
although the modelling approach used in that paper was com- 
pletely different. 

In our model we also assume a constant angle a for 
the direction of the velocity vector. In reality, this direction 
could vary with radius as the angular momentum distribution 
changes. For example, at smaller radii, a could be smaller as 
the velocity field more closely resembles pure Keplerian rota- 
tion. The inspiral time, and mass accretion rate, can therefore 
be very different. Only higher resolution observations can in- 
vestigate this further. 



Is L1489 IRS in any way special? No objects like it are re- 
ported in the literature. Since its life time is of the order of a 
few 10 4 yr, roughly 5-10% of the embedded phase, more ob- 
jects like it would be expected. It is possible that L1489 IRS 
was formed out of a core with unusually large angular momen- 
tum, which led to the formation of an untypically large disc. 
We cannot exclude that possibility, but hydrodynamical simu- 
lations are required to investigate how much angular momen- 
tum would be required, and if turbulent cloud cores can contain 
such amounts of rotation. 

An intriguing possibility is that the proximity of the neigh- 
bouring core is in some way related to L1489 IRS' special na- 
ture. The core is 60" (8400 AU) away in projection, and lo- 
cated in front of L1489 IRS, likely by a distance of comparable 
magnitude. Its systemic velocity is 0.4 km s _1 lower than that 
of L1489 IRS itself, indicating that L1489 IRS and the core 
currently are moving away from one another, at least in the di- 
rection along the line of sight. Could the neighbouring core be 
feeding material onto the disc of 11489 IRS? The velocity gra- 
dient is such that it merges smoothly with the velocities in the 
core (see Fig.0. This suggest that a physical link between the 
core and the disc may exist. If the core feeds material onto the 
disc, this could be significant source of angular momentum, 
thus keeping the disc large. It is, however, not easy to under- 
stand why gas in the core would be gravitationally bound to 
the L1489 IRS star, because of the significant mass reservoir in 
the core itself of 2.9 M . 

Another hypothesis is that L1489 IRS actually originated 
inside the core, but has since migrated away. Its current veloc- 
ity offset of 0.4 km s _1 is sufficient to move it to its current 
location 60" away in a few times 10 5 yr, the typical life time of 
an embedded YSO. We of course do not know how far its offset 
along the line of sight is, or what its three-dimensional velocity 
vector is like. Its line-of-sight velocity of 0.4 km s _! is not very 
different from the velocity dispersion of T Tauri stars and the 
turbulent motions in cloud complexes. One could propose that 
L1489 IRS' natal core was a turbulent, transient structure, and 
that, once formed, the YSO migrated with its gravitationally 
bound disc-like environment to a location outside the surround- 
ing cloud core, in effect 'stripping' the Class I object from most 
of its envelope. In this scenario we would now be seeing the in- 
ner, rotating Class I envelope around LI 489 IRS unobscured 
by the outer envelope. This might account for the different ap- 
pearance of L1489 IRS. 

This hypothesis can be tested in two ways. First, 'normal' 
Class I objects can be studied at high spatial resolution and 
in dense gas tracers to explore if the inner envelope is dom- 
inated by rotation on 1000 AU scales. Second, hydrodynam- 
ical calculations can be used to explore on what time scales 
newly formed YSOs can migrate away from their nascent cloud 
core. Searching for other objects like L1489 IRS would also 
be very useful. Potential targets would appear compact in sub- 
millimetre continuum images, but show strong emission lines 
in dense gas tracers. L1489 IRS shows HCO + lines with inten- 
sities of sever al K, while a fe w tenths of K is more typical for 
T Tauri discs dThi et alJ2004 . 
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5. Summary 

We have made a two-dimensional axi-symmetric disc-like 
model of the Young Stellar Object L1489 IRS. Line radia- 
tion transfer calculations produce synthetic spectra which can 
be directly compared to observations. We show that a flat- 
tened model gives a better description than a spherical one. 
Our model also reproduces millimeter interferometric imaging, 
near-infrared scattered light images, and CO ro-vibrational ab- 
sorption spectroscopy. 

We conclude further that the central star has a mass of 
1.35 M . The velocity vectors make an angle of 15° with the 
azimuthal direction. The velocity field is dominated by rota- 
tion, while small but significant amount of infall are present. 

A neighouring cloud core is present next to L1489 IRS. 
This cloud is well modelled by generic dark cloud parameters 
and we argue that it is likely situated in front of L1489 IRS. 
We speculate that this core could be feeding high angular- 
momentum material onto the L1489 IRS disc, explaining the 
unusual size of this object. Another explanation could be that 
L1489 IRS is an ordinary Class I object which has migrated 
away from its parental core, leaving it surrounded by only its 
gravitationally bound inner envelope. If this is true, L1489 IRS 
may provide valuable insight on the formation of protoplane- 
tary discs. 

This paper shows that it is possible to construct a global 
model of a Young Stellar Object that is able to fit observations 
on a wide range of spatial scales. Single dish line observations 
provide enough information to make highly detailed models 
of circumstellar structures, even on scales that are unresolved. 
However, on scales as small as 100 AU, the description may 
no longer be accurate, as evidenced by the near-infrared scat- 
tered light which suggest additional flattening of the disc-like 
structure. Recently obtained interferometer imaging with the 
Submillimeter Array may provide more insight on the inner 
several hundred AU around L1489 IRS (Brinch et al. in prep.). 
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